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Abstract 


This  paper  describes  a  series  of  numerical  simulations  of  field 
reversal  in  the  Reversed  Field  Pinch,  using  an  incompressible  MHD  model 
and  a  reference  set  of  plasma  conditions.  Field  reversal  and  maintenance 
are  observed,  although  at  values  of  the  pinch  parameter  8  larger  than  in 
experiment.  This  discrepancy  is  shown  to  arise  largely  from  the 
unrealistic  resistivity  profile  in  the  reference  conditions  and  may  not  be 
fundamental.  Qualitative  agreement  with  experiment  is  demonstrated  in 
several  areas.  The  view  that  field  reversal  is  due  to  a  simple  MHD  dynamo 


is  therefore  given  support. 
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Introduction 


This  paper  describes  the  results  from  a  series  of  numerical 
simulations  of  field  reversal  in  the  Reversed  Field  Pinch  (RFP),  using  a 
cylindrical  incompressible  magnetohydrodynamic  (MHD)  model.  The 
calculations  presented  here,  the  first  in  a  program  of  numerical  pinch 
modelling,  apply  specifically  to  the  plasma  conditions  chosen  in  [1]  as  a 
basis  for  numerical  comparison  work. 

The  purpose  of  the  paper  is  to  record  typical  simulation  results  for 
simple  macroscopic  quantities  of  experimental  interest  and  thus  help 
assess  to  what  degree  single-fluid  theory  can  explain  the  remarkable 
phenomenon  of  field  reversal.  These  calculations  extend  those  at  present 
in  the  literature  by  providing  a  sequence  of  results  rather  than  single 
instances  of  reversal. 

The  main  feature  of  the  present  MHD  model  is  the  assumption  of 
incompressibility.  This  was  adopted  as  a  physics  choice  as  a  particularly 
simple  model  with  which  to  explore  the  hypothesis  that  field  reversal  is  a 
gross  electromechanical  effect,  not  essentially  altered  by 
compressibility,  transport  etc. 

The  assumption  of  incompressibility  also  simplifies  the  numerical 
problem  since  magnetosonic  waves  and  their  severe  timestep  restrictions 
are  absent  from  the  system. 
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It  should  be  noted  that  the  plasma  conditions  taken  from  [  1  ] ,  notably 
the  low  value  of  Lundquist  number  S  and  the  flat  resistivity  profile, 
are  somewhat  different  from  experiment.  Furthermore,  numerical  simulation 
results  are  found  to  be  fairly  sensitive  to  the  choice  of  resistivity 
profile.  The  results  of  this  paper  should  therefore  be  clearly  understood 
as  applying  to  the  numerical  reference  simulation  proposed  in  [1],  rather 
than  an  attempt  to  model  precise  experimental  conditions. 

With  experiments  in  mind,  the  most  obvious  feature  of  the  present 
results  is  the  rather  large  value  of  pinch  parameter  6  required  for 
reversal  (6  *  2.3).  The  result  in  [4]  and  preliminary  results  with  the 
present  code  indicate  that  (even  with  an  incompressible  MHD  model)  at 
larger  values  of  S  and  with  a  resistivity  profile  increasing  towards  the 
edge  (nearer  real  conditions),  the  value  of  6  required  for  reversal 
drops  to  1.4  -  1.6,  much  more  consistent  with  experiment. 

Background 


Simulations  of  field  reversal  are  not  new  and  various  simulation 
codes  have  been  used  during  the  development  of  the  subject. 

The  first  3-dimensional  numerical  simulation  of  the  pinch  was  that  of 
Sykes  and  Wesson  [2] ,  using  a  compressible  MHD  model  with  anisotropic 
resistivity  in  cartesian  geometry.  This  demonstrated  field  reversal  and 
sustainment  arising  from  gross  MHD  activity,  but  did  not  identify  the 
precise  dynamical  processes  responsible. 
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Park  and  Chu  (3}  used  a  2-dimensional  incompressible  model  with 
anomalous  resistivity.  Field  reversal  was  obtained  and  a  successful 
comparison  made  with  the  field  profile  in  the  eta-beta  experiment. 

Aydemir  and  Barnes  [4]  used  a  3-dimensional  incompressible  MHD  model 
with  a  resistivity  profile  increasing  sharply  at  the  plasma  edge  and 
obtained  a  steady  reversed  state  with  8  =  1.5,  F  «  -0.03  at  S  «  5.103. 

Sato  [5]  has  proposed  a  driven  reconnection  model  of  field  reversal 
and  has  performed  calculations  with  a  compressible  MHO  model. 

In  [1],  a  comparison  of  field  reversal  calculations  was  undertaken 
using  compressible  and  incompressible  3D  MHD  codes.  This  clearly 
demonstrated  a  difference  between  the  compressible  and  incompressible 
results.  For  the  single  specific  case  run,  the  compressible  codes  gave 
field  reversal  whereas  the  incompressible  codes  did  not.  The  macroscopic 
reason  [1]  for  this  difference  appears  to  be  a  significant  contribution  of 
v  B  to  the  (v  *  B).  term  driving  the  j.  current  responsible  for 
field  reversal.  (Here  the  bar  denotes  the  axisymmetric  part.)  This 
contribution  is  identically  aero  in  an  incompressible  fluid  since  vf  =  0. 
However,  the  contribution  must  also  vanish  at  the  reversal  point  (5z  ”  0) 
and  therefore  cannot  be  solely  responsible  for  reversal,  irrespective  of 
the  plasma  mjdel.  The  significance  of  this  contribution  and  hence  the 
degree  of  applicability  of  an  incompressible  MHD  model  awaits  a  more 
detailed  comparison  with  experiment  and  a  more  detailed  study  of  the 
results  from  compressible  codes.  There  is  no  question,  however,  that 
field  reversal  can  be  obtained  with  an  incompressible  MHD  model. 


Numerical  model 


The  numerical  model  [6]  is  that  of  single  fluid  incompressible  MHO  in 
a  periodic  cylinder.  In  dimensionless  form  the  equations  may  be  written 

dv 

_  +  v.  V  v  -  -  Vp  +  j  *  B  +  v  V2  v 

dt 

SB 

—  -  Vx(vxB  -  _I!j) 

dt  S 

V.v  «  0 

h  -  T)(r) 

v  «  const. 

Here,  times  are  normalized  to  the  Alfv£n  time  based  on  the  (minor)  radius 
of  the  cylinder.  S  is  the  ratio  of  resistive  diffusion  time  to  Alfv£n 
time  and  ti  is  the  resistivity  profile.  For  the  present  case  S  «  10 3 
and  t)  «  1  [11.  The  length  of  the  cylinder  is  2it  (unit  aspect  ratio). 

Inclusion  of  the  viscosity  term  is  necessary  to  ensure  non-linear 
numerical  stability  of  the  code.  For  the  present  case,  it  is  sufficient 
to  take  v  »  10”1*. 

These  equations  are  advanced  in  time  using  an  explicit  2-step  scheme 
with  an  implicit  treatment  of  resistivity.  A  spectral  method  is  used  in 
the  6  and  z  directions  and  conventional  second  order  finite 
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differences  in  radius.  Incompressibility  is  enforced  by  use  of  the 
standard  MAC  method.  A  detailed  description  of  the  code  is  available. 

The  boundary  conditions  are  those  pertaining  to  an  impermeable, 

perfectly  conducting,  free-slip  wall  with  a  constant  uniform  driving 

electric  field,  i.e.  v  »  0,  B  -  0,  u.  ~  0,  u  •  2v7r,  E.  «  0, 
r  r  e  z  er  e 

E  «  const.  (Here  oj  =  V  x  v.) 
z  -  - 

The  initial  magnetic  field  is  taken  to  be  the  force-free  equilibrium 
with  safety  factor  profile  q  given  by  [1] 

q(r)  -  0.4  (1.0  -  1.8748  r2  +  0.83232  r1*). 

(The  field  profiles  are  then  obtained  numerically.)  This  means  that  the 

initial  conditions  have  a  reversed  B  .  It  must  be  stressed  that  the 

z 

simulations  of  field  reversal  described  here  are  not  a  consequence  of  the 
choice  of  a  reversed  initial  state.  Field  reversal  is  certainly 
obtainable  starting  from  the  more  realistic  case  B^-  conBt.,  which  indeed 
is  the  initial  condition  normally  used  in  the  code.  The  initial  velocity 
field  is  zero  everywhere. 

The  equilibrium  magnetic  field  is  seeded  with  white  noise  at  a 

relative  level  of  10-6  and  the  code  run  with  E„  «  0  (constant 

'toroidal*  flux)  and  E  -  const,  for  a  substantial  fraction  of  a  field 
z 

diffusion  time.  The  value  of  E  is  chosen  to  give  an  appropriate  value 

z 

of  the  pinch  parameter  6  in  the  final  state  (and  may  therefore  be 
different  from  that  implied  by  the  initial  q  profile). 
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The  results  described  here  were  obtained  using  a  modest  computational 

mesh  of  N  -  32,  N.  -  12,  N  “  25,  This  means  that  pololdal  mode  numbers 
r  9  z 

0  <  m  <  3  and  toroidal  mode  numbers  -8  <  n  <  8  are  Included  In  the 
calculations. 

Results 


Fig.  1  shows  a  typical  plot  of  the  pinch  parameter  8  and  field 
reversal  fraction  F  as  a  function  of  time.  As  the  simulation  proceeds 
the  field  rapidly  loses  the  reversal  of  the  initial  state,  but  then 
suddenly  reverses  and  remains  marginally  reversed  for  the  remainder  of  the 
simulation,  although  with  slight  fluctuations.  This  remarkable  behaviour 
is  emphasized  by  the  dashed  lines  showing  how  F  and  6  would  evolve 
(under  axisymmetric  resistive  diffusion)  if  there  were  no  dynamo  effect. 

Fig.  2  shows  F  versus  0  for  the  present  plasma  model  (labelled 

S  «  10 3  T)  »  flat).  It  is  obtained  from  a  sequence  of  separate  simulations 

(each  like  Fig.  1)  with  different  values  of  the  applied  electric  field 

E  .  Here  and  below  we  use  the  convention  that  each  point  on  the  graph  is 
z 

the  time-averaged  value  in  the  respective  quasi  steady-state.  The  error 
bars  then  indicate  the  maximum  deviation  over  time  from  this  average.  The 
diagram  clearly  demonstrates  that  field  reversal  can  be  obtained  with  an 
incompressible  MHD  model.  However,  for  the  present  case,  the  value  of 
pinch  parameter  8  required  for  reversal  (8  =  2.3)  is  larger  than  that 
observed  in  experiment  (8  =  1.4).  This  result  and  its  known  dependence  on 
the  choice  of  plasma  conditions  have  been  discussed  above. 
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To  illustrate  this  point,  Fig.  2  also  shows  results  obtained  from  a 
more  realistic  plasma  model  (labelled  S  »  1 0 4  t)  =  shaped).  This  has  a 
resistivity  profile  increasing  towards  the  edge  (the  most  significant 
change)  and  a  larger  S  value.  For  this  case,  the  F  -  8  curve  lies 
considerably  to  the  left  of  that  for  the  present  case  of  flat  resistivity 
and  is  much  more  consistent  with  experiment.  It  is  also  consistent  with 
the  result  in  [4] . 

Fig.  3  shows  the  corresponding  r.m.s.  values  of  the  relative 

fluctuation  level  6B  /B  at  the  wall  plotted  against  the  driving 
0  0 

electric  field  E  .  (Here  6B  =  B  -  B  .)  The  typical  fluctuation 
z  0  0  0 

level  of  ^  io%,  rather  larger  than  the  experimental  value  of  =  1%,  is  also 
likely  to  be  dependent  on  the  choice  of  plasma  conditions.  Nevertheless, 
it  is  significant  that  the  fluctuation  level  is  not  100%,  as  might  be 
expected  from  simplistic  quasi-linear  MHD  scaling  arguments.  Whether  MHD 
models  always  predict  a  relatively  large  level  of  fluctuation  and  thus 
whether  a  more  complicated  description  is  needed  in  this  respect  [7]  can 
be  decided  only  by  further  simulation  work. 

Fig.  4  shows  the  resistance  per  unit  length  (E  /I)  as  a  function  of 

z 

.  Shown  for  reference  is  the  corresponding  resistance  obtained  if  the 
plasma  acted  as  a  solid  conductor  of  the  same  resistivity.  The  clear 
non-linear  relationship  between  the  current  and  the  applied  voltage 
constitutes  a  numerical  simulation  of  what  is  known  experimentally  as 
current  'screw-up'. 


Fig.  5  shows  simulation  results  for  the  p  profile  (the  ratio 


j  /B)  as  a  function  of  radius  and  makes  a  heuristic  comparison  with 

a 

experimental  values  from  the  HBTX  1A  experiment  [8] .  The  simulation  curve 
corresponds  to  the  case  of  greatest  reversal  in  Fig.  2  (S  =  10  3,  T)  ■=  flat, 
8  “2.9,  F  =  -0.21)  and  was  calculated  at  the  end  of  the  run  (time 
t  «  320).  It  has  a  value  of  F  nearest  that  of  the  experimental  case. 

The  experimental  curve  was  obtained  for  S  >  10**,  6  “  1.7,  F  =  -0.3, 
aspect  ratio  =  3  and  is  the  only  direct  measurement  of  p  available  from 
HBTX  1A.  Although  the  differences  in  the  cases  preclude  a  true 
comparison,  there  is  an  encouraging  consistency  between  the  simulation  and 
experimental  results. 

Fig.  6  shows  the  radial  profile  of  the  axisymmetric  part  of  j  at 

b 

the  end  of  the  calculation  for  the  case  8  “  2.3,  F  “  0.  Of  fundamental 
physics  interest  is  the  macroscopic  cause  of  this  current,  responsible  for 
the  decreasing  profile  (through  jQ  “  -  dB^/dr).  To  address  this 

point,  the  axisymmetric  parts  of  Ohm's  law 

—  j„  “  E„  +  <(v  x  B)  „> 

e  8  -  -  9 

are  also  displayed  on  the  diagram.  It  is  clear  that  <vzBr>  is  the 
principal  driving  term  for  the  field  reversal,  <vrBz>  is  a  minor 
contribution  and  Eg  =0.  The  importance  of  <vzBr>  has  also  been 
suggested  by  Sato  [5] • 

Fig.  7  shows  a  logi Q  plot  of  the  energy  spectrum  (the  energy  in  each 
Fourier  harmonic)  of  B^  as  a  function  of  time  for  the  case  6  “2.3, 

F  «  0.  (For  reference,  the  corresponding  log10  energy  spectrum  value  of 


Bg  is  typically  +0.1.)  It  shows  the  dominance  of  the  two  inodes  m  =  1, 
n  “  -2  and  m  =  1,  n  =  -3  throughout  the  simulation.  (For  simplicity  of 
representation,  only  these  modes  have  been  labelled  on  the  figure.)  In 
general,  the  dominant  modes  are  found  to  be  resonant  inside  the  reversed 
surface.  The  importance  of  m  =  1  modes  is  consistent  with  experiment 
[9)  • 


Fig*  8  shows  the  energy  spectrum  at  the  end  of  the  same  calculation 
(time  t  *  500)  as  a  histogram  in  (m,n)  space.  It  provides  a  simple 
picture  of  the  amplitude  of  the  various  modes,  in  particular  the 
non-linear  m/n  sequence  1/-2,  1/-3,  2/-5,  3/-8. 

Concluding  Remarks 


3-dimensional  numerical  simulations  of  field  reversal  in  the  Reversed 
Field  pinch  have  been  performed  using  (as  a  physics  choice)  the  simplest 
incompressible  MHD  model  and  a  reference  set  of  plasma  conditions  [1). 
Field  reversal  and  maintenance  are  clearly  observed.  The  F  -  9  curve 
has  been  presented. 

There  is  considerable  qualitative  agreement  between  the  numerical 
results  (F  -  6  curve,  fluctuation  levels,  fluctuation  spectrum,  current 
screw-up  and  p  profile)  and  experimental  observations. 

For  the  present  case,  the  numerical  values  of  the  pinch  parameter  8 
and  fluctuation  levels  are  rather  larger  than  in  experiment.  However, 
this  quantitative  discrepancy  is  lenown  to  arise  largely  from  the  choice  of 
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somewhat  unrealistic  plasma  conditions,  and  may  simply  reflect  the  need 
for  correspondingly  increased  MHD  activity. 

Taken  together,  these  results  lend  weight  to  the  view  that  field 
reversal  is  a  phenomenon  arising  principally  from  single  fluid  MHD  dynamo 
action. 

In  common  with  all  previous  work,  the  present  simulations  are 
essentially  a  series  of  macroscopic  numerical  experiments  and  do  not 
attempt  to  identify  the  detailed  nature  of  the  modes  responsible  for  the 
dynamo.  Calculations  are  in  hand  to  explore  this  question,  to  establish 
the  results  for  more  realistic  S  values,  aspect  ratios  and  resistivity 
profiles,  and  to  make  further  comparisons  with  experiment. 
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F  and  6  vs.  time.  Dashed  lines  show  evolution  under  ID  resistive 
diffusion.  (Case:  S  *  10 3,  t\  -  flat). 


0  12  3 

F  vs.  0  in  the  steady  state.  Bars  show  fluctuation  amplitudes 
(Cases:  S  -  10 3,  t|  -  flat  and  S  -  10 4,  t)  -  shaped). 


6B g/Bg  vs.  Ez  in  the  steady  state.  Bars  show  fluctuation 
amplitudes. 


(case:  S  «  10 3,  ti  -  flat). 


Fig. 4  vs*  Er  in  the  steady  state.  Bars  show  fluctuation 


amplitudes. 

(Case:  S  -  10 3,  r)  -  flat). 
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Fig. 5  j(  /B  vs.  radius  and  hauristic  comparison  with  axpsrimai 
(Case:  S  «  103,  r)  -  flat). 


Fig. 6  j  and  (S/r))  (  v  x  B)  terms  vs.  radius. 
0  “*  0 

(Case:  S  -  103,  f)  ■  flat). 
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